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O ' ABSTRACT 

<n : 

We investigate the nonlinear two- and three-point correlation functions of the cosmo- 
' logical density field in Fourier space and test the popular hierarchical clustering model, 

that the three-point amplitude Q is independent of scale. In high-resolution TV-body 
y—i ■ simulations of both n = —2 scale-free and cold dark matter models, we find that Q at 

late times is not constant but increases with wavenumber far into the nonlinear regime. 
Self-similar scaling also does not hold as rigorously for the three-point function as for 
the two-point function in the n = —2 simulation, perhaps a manifestation of the finite 
simulation volume. We suggest that a better understanding of the two- and three-point 
correlation functions in the nonlinear regime lies in the link to the density profiles of 
' dark matter halos. We demonstrate and quantify how the slopes of the correlation func- 

tions are affected by the slope of the halo profiles using simple halo shapes and analytic 



clustering models. 



Subject headings: cosmology : theory — dark matter — large-scale structure of universe 
— methods: numerical 



1. Introduction 

In this Letter we examine the power spectrum and the bispectrum for numerical simulations 
of the n = —2 scale free model and cold dark matter (CDM) model. The power spectrum P(k) 
is the most fundamental statistic for a random field and the only statistic for a Gaussian field. 
The cosmological distribution of matter, however, is not Gaussian, and higher order statistics 
are also important. The lowest-order non-Gaussian information is encoded in the bispectrum 
B(ki,k2,ks). In analogy to P(k), which is the Fourier transform of the two-point correlation 
function £ , the bispectrum is the Fourier transform of the three-point correlation function £ . The 
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power spectrum P{k) of matter fluctuations is related to the mass density field 5 = p/p — 1 in k- 
space by (5(ki)5(k 2 )) = P(ki) 8r>(ki + k 2 ) > where 5r> is the Dirac delta-function. The bispectrum 
B(k\, k 2 , k 3 ), similarly, is related to 5 by (5(k 1 )5(k 2 )5(k 3 )) = B(ki,k 2 , fc 3 ) <5d(^i + k 2 + k 3 ) . 

On sufficiently large length scales, density fluctuations are small enough that perturbation 
theory for gravitational instability suffices, and the predictions for the correlation functions can 
be computed in detail. To linear order, the fluctuation amplitude 5(k,t) grows by an overall scale 
factor, 5 = D(t)5o(k), where 6o(k) is the amplitude at some early time to and D(t) is a growing 
function of time. In a nonlinear theory, interactions generate a sum of terms to all orders in 
the initial amplitude, 5 = S^ 1 ' + 5^ + 5^ + • • • , where 5^ ~ 5 n . Even for a Gaussian initial 
distribution, these nonlinear terms induce nonvanishing higher order correlations. In perturbation 
theory, the lowest order two-point statistic is the linear power spectrum Pi(k) oc [t^ 1 - 1 ] 2 . For the 
three-point statistic, the first nonvanishing contribution (i.e. tree-level) to the bispectrum is related 
to Pi by 

B (0) (fci, h, h) = F(k u k 2 ) Pi{k{)Pi{k 2 ) + F(k 2 , k 3 ) PiihWh) + F(k 3 , fci) Pik^Pih) , (1) 

where F(k u kj) = ^ + (h/kj + kj/h) {hi ■ kj) + f (fei • kj) 2 (Fry 1984). For Q < 1 the factors ^ 
and I become 1 ± k, where for Q near 1, k ~ | in an open universe and k ~ | O -1 / 143 in a 

flat universe with cosmological constant (Bouchet et al. 1995). 

Equation (|l]) indicates that the main dependence on the power spectrum can be removed if 
one defines the hierarchical three-point amplitude 

Q[k u fc 2; fc 3 j _ p(A;i)iW + P(A;2) p (A;3) + p WP(il) • 

The three-point amplitude Q has the convenient feature that for the lowest nonvanishing result in 
perturbation theory (eq. Q), is independent of time and the overall amplitude of Pi. Moreover, 
for scale- free models with a power-law Pi(k), is independent of overall scale as well. In a pure 
hierarchical model, Q is exactly constant. 

The bispectrum B and the three-point amplitude Q depend on any three parameters that define 
a triangle in /c-space. For simplicity, we will use equilateral triangles to explore scale dependence. 
We have examined several other configurations and found behavior similar to the equilateral case. 
For a closed equilateral triangle, one has k\ = k 2 = k 3 = k and ki ■ kj = — |, and the bispectrum 
B eq depends only on the single wavenumber k. To lowest order, the equilateral bispectrum has a 
particularly simple form, Beq (k) = ^-Pf(k) , and it follows from equation (]2|) that (^^{k) = | , 
independent of the power spectrum. Since B has dimensions of Mpc 6 , it is convenient to define 
a new dimensionless variable, which we name A3, to characterize the three-point statistic of mass 
density fluctuations. Following the familiar density variance A (A;) = 47rfe 3 P(k) for the two-point 
statistic, we define A3 for equilateral triangles as 



A 3 (fc) = 4vrA: 3 J B cq (k) ; (3) 
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and it follows that As(k)/A(k) = ^3 Q cq (k). Since Q cq is simply a constant to the lowest nonva- 
nishing order in perturbation theory, this new statistic also has the convenient property that the 
leading order Af\k) and A(°)(ife) have identical shapes, i.e., Af\k)/A^(k) = (12/7) 1 / 2 . 



2. High-/c Behavior 

The leading-order perturbative expressions given above are valid for i < 1, or at low k. 
Including the next-to-leading order terms for the bispectrum extends the range of validity to the 
quasilinear regime 8 ~ 1 (Scoccimarro et al. 1998). Numerical simulations, however, are required 
in order to determine the fully evolved bispectrum of mass density fluctuations into the deeply 
nonlinear regime at high k. 

We examine results from two cosmological simulations, an n = — 2 scale-free model and a 
cold dark matter (CDM) model. The n = — 2 simulation has 256 3 particles and a Plummer force 
softening length of L/5120, where L is the box length. This is the highest resolution simulation of 
an n = —2 model to our knowledge. It is also the main simulation used in the study of self-similar 
scaling of the power spectrum by Jain & Bertschinger (1998). The CDM model is critically flat 
with matter density Q m = 0.3 and cosmological constant S7a = 0.7. This run has 128 3 particles and 
is performed in a (100 Mpc) 3 comoving box with a comoving force softening length of 50kpc for 
Hubble parameter h = 0.75. The baryon fraction is set to zero for simplicity. The primordial power 
spectrum has a spectral index of n = 1, and the density fluctuations are drawn from a random 
Gaussian distribution. The gravitational forces are computed with a particle-particle particle- 
mesh (P 3 M) code (Ferrell & Bertschinger 1994). We compute the density field 5 on a grid from 
particle positions using the second-order triangular-shaped cloud (TSC) interpolation scheme. A 
fast Fourier transform is then used to obtain 5 in fc-space. The /c-space TSC window function is 
deconvolved to correct for smearing in real space due to the interpolation, and shot noise terms are 
subtracted to correct for discreteness effects. 

Figure 1 illustrates the behavior of the nonlinear power spectrum A(k) and the three-point 
statistics A^(k) and Q(k) for equilateral triangular configurations for the n = —2 scale-free simu- 
lation. Three time outputs are shown, and the results are plotted against scaled k/k n i, where k n i 
is defined by d 3 k Pi(a,k) = 1, and the expansion factor a is 1 initially. Unlike the two-point 
A, for which self-similar scaling works well (Jain & Bertschinger 1998), the three-point curves in 
the plot do not overlap; self-similarity therefore does not hold as rigorously for the three-point 
function in the simulations. The slope of Q also has a weak dependence on time. For the last 
output (a = 26.91, k n i = 7.25), Q increases monotonically with k down to the resolution limit with 
an approximate logarithmic slope of 0.225. The next-to-leading order perturbative results for Q 
(Scoccimarro et al. 1998) roughly agree with the most weakly evolved, (a = 13.45, k n i = 29) output 
up to k ~ k n i, but not at higher k nor in the later outputs. The upturn of Q cq at high k (dotted 
portion of curves in Fig. 1) is related to the downturn in A at the same k, which is likely due to 
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limited numerical resolution. 

Many earlier papers did not report this behavior (e.g., Efstathiou et al. 1988; Pry, Melott, 
& Shandarin 1993; Scoccimarro et al. 1998). Some of the computations were not carried out 
to sufficiently high wavenumbers to provide a clear test of the hierarchical model. Others (e.g. 
Scoccimarro et al. 1998) assumed self-similarity and averaged results over different time outputs, 
treating each output as if it were from an independent realization. We have checked for possible 
numerical artifacts by re-examining the n = — 2 results of Fry et al. (1993) from 128 3 particle-mesh 
(PM) simulations. At evolution stages comparable to those in Figure 1 (i.e. k n \ = 32, 16, and 8), 
we find that Q eq in four different random realizations agrees well with Figure 1 up to the resolution 
limit of the PM runs. 

Departures from the hierarchical model have been suggested in several papers that have ex- 
amined the real-space three-point function (Suto & Matsubara 1994) and the skewness S3 (e.g., 
Lahav et al. 1993; Colombi et al. 1996), a volume integral average of the full three-point function. 
The simulations used in these papers had 64 3 particles and a smaller dynamic range, so the results 
were limited to lower k. The skewness is also a less powerful probe of hierarchical clustering than 
the full three-point examined here because the volume integral might average out non-hierarchical 
signals. 



3. Dependence on Halo Density Profiles 

In this section we examine analytical models that can be used to elucidate the numerical 
results of Figure 1. We suggest that the key to a better understanding of Q in the high- A; regime 
lies in the connection between Q and the density profiles of dark matter halos. This is because the 
contributions to Q at high k are dominated by close particle triplets, many of which are located 
within individual halos of size ~ 27r/k. It is therefore reasonable to expect the two- and three-point 
correlation functions on small length scales to reflect the density profiles of dark matter halos. 
Sheth and Jain (1997) have examined how the amplitude of the two-point function is related to 
simple halo profiles. Our emphasis here is on the link between the shape of the three-point function 
and halo profiles. 

One simple way to demonstrate the relation between halos and the correlation functions is to 
consider the "power law cluster" model. In this model, matter is distributed in randomly placed 
halos with a power-law mass density profile p{r) oc r~ e out to some radius R (Peebles 1974; 1980). 
For 1.5 < e < 3 and r < R, one can derive analytically that the two-point function £, the density 
variance A, the three-point function (, and Q scale as 

f(r) oc r 3 - 2e , A(fc) oc k 2e ~ 3 , £(r) oc r 3 " 3e , Q oc ^ oc r e ~ 3 oc k 3 ~ e . (4) 

This model is clearly simplistic because in realistic cosmological models, the centers of dark matter 
halos are not randomly distributed, the halo density profile is not a pure power law, and not 
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all mass is located within the halos. Equation (Q) nonetheless provides valuable insight into the 
link between the shape of dark matter halos and the high-A; behavior of the two- and three-point 
clustering statistics. 

To test the applicability of equation (Q) in a more realistic setting, we have performed a series 
of experiments with dark matter halos identified in cosmological iV-body simulations. In these 
experiments, we keep the centers and masses of the halos unchanged but redistribute the subset of 
particles which lies within the virial radius (the radius within which S > 200) of each halo according 
to p oc r~ e for a chosen e. We then recompute A and Q eq from the particle positions. This recipe 
allows us to explore simple density profiles while preserving the actual halo-halo spatial correlations, 
the halo mass distribution, and the correlation signal from non-halo particles in the simulations. 
Figure 2 compares the results for the original output in the n = — 2 scale- free simulation with the 
results for replaced power-law halos. Figure 3 does the same for the Q m = 0.3 CDM simulation. 
As indicated, when the halo profiles obey p oc r _e , we find that A and Q at high k are well 
approximated by power laws, and that the slopes are indeed related to e according to equation (Q). 
For e = 2 (i.e. an isothermal density distribution), for example, both Figures 2 and 3 show that 
A oc k and Q eq oc k. Figure 3 also verifies equation (Q) for e = 2.35: A oc k 1,7 and Q cq oc /c 65 . 
We therefore conclude that for a power-law density profile, equation (Q) works well at predicting 
the slopes of the power spectrum and the three-point Q in the nonlinear regime. This remains true 
even when the clustering of the halo centers and the correlation signal from non-halo particles are 
taken into account. 

4. Discussion 

To gain a better understanding of the behavior of two- and three-point statistics in the deeply 
nonlinear regime, we have investigated in this Letter how the shapes of dark matter halos affect the 
power spectrum P and the three-point amplitude Q at high k. For halos with power-law density 
profiles, Figures 2 and 3 demonstrate that the simple model of equation (||) works very well at 
predicting the slopes of P and Q at high k. The figures also show that the converse is not true: the 
slope of P or Q alone does not determine a unique halo shape. For example, Figure 2 shows that P 
at high k in the n = — 2 scale free model is only slightly changed when simulation halos are replaced 
with e = 2 power-law halos. It would therefore be difficult to distinguish these two profiles using 
the power spectrum alone. The three-point Q, however, is drastically different for the two profiles. 
This indicates that although the nonlinear P and Q are both closely related to the shapes of halos, 
they probe different regions and properties of the halos. A plausible explanation for this difference 
is that the two-point function reflects the mean density at a given distance from a particle, whereas 
the three-point function reflects the dispersion in density at the given scale (Peebles 1980, § 36). 

Realistic dark matter halos do not have pure power law profiles. High-resolution iV-body 
simulations have shown that cold dark matter halos have a roughly (Jing & Suto 2000) universal 
density profile with e = 3 at the outer radii and e = 1 (Navarro, Frenk & White 1997; Huss, Jain, & 
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Steinmetz 1999) or e = 1.5 (Moore et al. 1999) nearer the halo centers. A further complication is 
that these profiles have a mass-dependent "concentration parameter" : less massive halos are more 
centrally concentrated. It will be useful to develop a model more extensive than equation (Q) to 
incorporate these factors. Such a model should help to shed light on why Q(k) at high k in our 
figures is approximately a power law, even though the halo profile in the simulations is not a pure 
power law. Given that the shape of Q is highly sensitive to halo profiles, the model can also be 
used to tackle the question: what criteria must the halos satisfy in order for hierarchical clustering 
to occur? 

In the best simulation of the n = —2 scale-free model available to us (with 16.7 million 
particles and a force resolution of L/5120), we find that the three-point amplitude Q is not a 
constant with wavenumber k and that the self-similar scaling observed for the two-point function 
(Jain & Bertshinger 1998) does not hold as rigorously for the three-point function (Fig. 1). We 
have checked for possible numerical artifacts using smaller simulations with different phases (see 
§ 2), and the results are similar within the range of overlap but may still be a manifestation of 
the finite simulation volume, to which the n = —2 spectrum is particularly susceptible. Based on 
lower resolution simulations for scale-free initial conditions with spectral index n, Fry et al. (1993) 
suggested Q = 3/(3 + n) as the asymptotic value in the nonlinear regime. Scoccimarro & Frieman 
(1999) in so-called hyperextended perturbation theory obtain Q = (4 — 2 n )/(l + 2 n+1 ). These two 
expressions would predict Q = 3 and 2.5, respectively, for n = — 2 at high k. Figure 1 shows that 
Q reaches a similar value between 2.5 and 3 at k ~ k n i, but it rises to 4 to 5 at k ~ 10 k n i. For 
CDM-type models, the effective index decreases with k, so the two expressions above would give Q 
rising with k, but in neither case is the prediction as strong as what we see. 

It is useful to compare our numerical results with theoretical expectations. From the as- 
sumptions of stable clustering and self-similarity, it is well known that for the two-point function, 
A(k) oc fc(9+3n)/(5+n) a t high k for an initial P t oc k n (Davis & Peebles 1977). The nonlinear A in 
Figure 1 and earlier work of Efstathiou et al. (1988), Peacock & Dodds (1996), Jain (1997), Ma 
(1998), and Jain & Bertschinger (1998) all support this behavior. For the three-point function, 
these same conditions would result in constant Q, i.e., hierarchical clustering. The fact that high- 
resolution simulation results do not support this implies that one of these assumptions - stability 
or self-similarity - must be violated, or still better future simulations are needed. 

We are grateful to John Peacock for enlightening discussions and Edmund Bertschinger for 
providing the n = —2 scale-free simulation. Computing time for this work is provided by the 
National Scalable Cluster Project and the Intel Eniac2000 Project at the University of Pennsylvania. 
C.-P. M. acknowledges support of an Alfred P. Sloan Foundation Fellowship, a Cottrell Scholars 
Award from the Research Corporation, a Penn Research Foundation Award, and NSF grant AST 
9973461. 
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Fig. 1. — Scaling of the power spectrum, A(fc) = 47r/c 3 P(/c), and the three-point statistics A^(k) = 
Airk 3 Bll 2 (k) and Q C q(k) for equilateral triangles, in an n = —2 scale-free simulation. Three 
time outputs are shown, where the expansion factor (1 initially) and nonlinear wavenumber (in 
units of 2vr/L) are: (a,k n i) = (13.45, 29), (19.03, 14.5), and (26.91,7.25) (left to right). In the top 
panel, the solid and dashed curves are A and A3, respectively, computed from the simulation; the 
dotted straight line is the linear Ai ex k; and the dotted curve is the analytic approximation for 
the nonlinear A from Peacock & Dodds (1996). In the bottom panel, the dotted curves indicate 
the leading (horizontal line) and next-to-leading perturbation theory results for Q eq . The dotted 
portions of the curves at high k indicate regions of possible spurious effects due to finite numerical 
resolution. 
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Fig. 2. — The power spectrum A(k) = 4-7r/c 3 P(/c) and the three-point amplitude Q C q{k) for the last 
output {k n i = 7.25) of the n = — 2 model in Fig. 1. In both panels, the solid curves are computed 
from the original simulation particles. The dashed curve compares the results when the simulation 
particles in halos are redistributed according to pure power-law p oc r~ € with e = 2. The straight 
dotted line marks the slope predicted by eq. (Q), which works remarkably well at k ;> 70. 
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Fig. 3. — Similar to Fig. 2, but for a low-density CDM model with O m = 0.3 and Q\ = 0.7. In 
both panels, the solid curves are computed from the original simulation particles. The dotted curve 
tracing the solid one in the top panel shows the analytic approximation for A from Ma (1998). The 
other curves compare the results when the simulation particles in halos are redistributed according 
to pure power-law p oc r~ € with e = 2 (long-dashed) and e = 2.35 (short-dashed). The straight 
dotted lines mark the slopes predicted by eq. (||), which works remarkably well at k ^ 3/iMpc -1 . 



